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ABSTRACT 

We assess the robustness of a low Mach number hydrodynamics algorithm for 
modeling helium shell convection on the surface of a white dwarf in the context 
of the sub-Chandrasekhar model for Type la supernovae. We use the low Mach 
number stellar hydrodynamics code, MAESTRO, to perform three-dimensional, 
spatially-adaptive simulations of convection leading up to the point of the ignition 
of a burning front. We show that the low Mach number hydrodynamics model 
provides a robust description of the system. 

Subject headings: convection - hydrodynamics - methods: numerical - nuclear 
reactions, nucleosynthesis, abundances - supernovae: general - white dwarfs 



Introduction 



Sub-Chandrasekhar models for Type la form an attractive progenitor candidate be- 



cause of the abundance o 
traces back to the work of 



low mass white dwarfs. The modern model of these explosions 



Livnd fll99nyLivne fc Clasned (1199(11 ): IWooslev fc Weaver! ( 119941 ): 



Wiggins fc Falld ( 119971 ): I Wiggins et al.l ( 119981 ) amongst others. In this "double detonation" 
model, a detonation ignited in the accreted helium layer on the surface of a low mass car- 
bon/oxygen white dwarf drives a shock inward, compressing the star and initiating a deto- 
nation in the core. At the time of these calculations, a common concern was that the lack of 
resolution in simulations prevented a realistic investigation of these events. Recently, how- 
ever, the sub-Chandrasekhar mass progenitor model of Type la supernovae has seen renewed 
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interest flFink et al.l 120071 l2010t ISim et al.ll2010t IShen et al.ll2010[ ). as observations show an 
increasing diversity of SNe la events. 

As in the traditional Chandrasekhar-mass model of SNe la, the ignition of a burning 
front is preceded by a long period of convection. Here, however, the convection takes place in 
the He layer, driven by reactions at the base of the layer. To date, numerical simulations of 
this type of explosion have been initialized by seeding a detonation in the He layer. However, 
it remains an open question as to whether the turbulent convective flow in the He layer can 
ignite a detonation in the first place; it is als o quite possible that th ere can be multiple 
nearly- simultaneous ignitions (see for example. iGarcia-Senz et al.lll999[ ). 



A potential issue with this model is that a detonation in the He shell would pro- 
duce large amounts of ^^Ni at the edge of the star that is inconsistent with observations 
f lHoefiich fc KhokhlovlEooa iHoefiich et allll996[ ) (also see iNugent et aUllDQTl for additional 
concerns about the spectra). The mass of the He layer is uncertain, and it has been sug- 
gested that one way to address the overproduction of Ni is for the He shell to b e very thin, 
in which case it may be able to detonate without over-producing surface ^^Ni (IFink et al. 



2OIOI : iKromer et al.ll2010l ). Detailed one-dimensional stellar evolution calculations suggest 
that the helium in these srn allest mass shells may never ignite as a detonation to begin 
with (jWoosley &: Kaserul201ll ). If it does ignite, iTownsley et al.l (120 12[ ) showed that a robust 
detonation can propagate through the layer. 

The dynamics of the convection and development of a burning front in the sub-Chandrasekhar 
model is inherently multidimensional. Because the fluid and flame velocities in this initial 
phase are much lower than the sound speed, we can apply the same low Mach number 



methodology used to study ignition in t he Chandrasekhar-mass scenario ( IZingale et al.ll2009 



Zingale et al.l 1201 ll : iNonaka et al.l 120121 ) . Once we ensure that we have a robust simulation 
methodology, there is a large parameter space of initial models to explore. The goal of this 
first paper is to demonstrate that low Mach number hydrodynamics provides an efficient and 
accurate simulation platform to explore the convective stage of these sub-Chandra events. 



2. Numerical Methodology 



To study the turbulent c onvection in the sub -Chandrasekhar model, we use the MAE- 
STRO code as documented in lNonaka et al.l (120101 ). For an overview of the low Mach number 



equat ions and numerical methodology, we refer the reader to Section 2.1 of iNonaka et al. 
(I2OI2I ). In summary, MAESTRO is a finite- volume, adaptive mesh stellar hydrodynamics 
code that models the fiow using a low Mach number approximation — sound waves are filtered 
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out of the system, but compressibility effects due to stratification and local heat release are 
retained. For most of the simulations presented here, the star is modeled on a Cartesian grid 
with the center of the white dwarf at the coordinate origin, i.e., we model one octant of the 
full star. Since the convection is confined to the outer shell, away from the center of the star, 
this geometry captures the convective behavior well. The base state pressure and density 
are represented by one- dimension al radial profi l es wh i ch capture the hydr o static state of 



the st ar. As in the simulations in IZingale et al.l (120091 ): IZingale et al.l ( l201ll ): iNonaka et al. 



( I2OI2I ). we derive the temperature from the equation of state given the pressure, density and 
mass fractions, rather than evolving the enthalpy equation. Below we provide additional 
details specific to the simulations presented here; these introduce a variety of parameters 
that will be explored in the results section. 



2.1. Microphysics 



We use a general, publicly available equation of state consisting of ions, ra diation, and ar 



bitrarily degenerate/rela tivistic electrons, together with Coulomb corrections ( lTimmesll2008 



Timmes fc Swestyll2000l ). We use a simple r eaction network cons i sting of the triple-alpha and 
^^C(a,7)^^0 reactions . The rates are frora ICaughlan fc Fowlerl (I19881). with screening as i n 
Graboske et~aD Jl973h : IWeaver et allJlQTsh - lAlastuev fc JancovicilJlQTshilltoh et al.Nl979h 
?he ^^ C(a, 7)^^0 reaction rate has been multiplied by 1.7 as suggested bv I Weaver fc Wooslev 
(I1993 ) and lGarnettI (119971 ). This network is an extension of the network used by lMalone et al 



(120111 ). When computing the effect of reactions over a time interval, we evolve the temper- 



ature along with the mass frac tions, keeping the the rmodynamic derivatives frozen during 
the integration, as described in lAlmgren et al.l ( l2008l ). 



2.2. Initial Model 

For simplicity, we construct our own one- dimensional, semi-analytic initial model of a 
white dwarf with a helium layer. This allows us to control specific features as we learn about 
the algorithmic sensitivity to the choice of initial model parameters. The initial model is 
constructed with an isothermal C white dwarf with an isentropic He layer on the surface. 
When initializing the data on the three-dimensional grid, we interpolate from this initial 
model and then apply a velocity perturbation, the latter described in Section 12.51 

To construct the initial model, we use the following iterative process: 
• We start by providing an estimate for the white dwarf central density, Pcore? and the 
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density at which we transition to hehum, puc- We also specify the composition of 
the core, Xcore, and its temperature, T^ore, and compute the central pressure via the 

equation of state, Pcorc = P(Pcorc, ^corc, -^corc)- 

We then construct our model by integrating outward from the center and iterate over 
the central density and transition density until the mass of the core, Mwd and the 
mass of the helium envelope, Mne, are the desired values. 

To be more precise, given pcore and pne, we do the following: 

— Specify the composition for zone i: 

{^core if Pi-1 > PHc 

Xcorc + K^He " ^core) [l + tauh (M^)] if i [l + tauh (^.^^H^)] < 0.999 
X}ie otherwise 

(1) 

here, 6 is the width of the transition layer and xhc = a;(pHc) is the coordinate 
corresponding to phc- It is important that 6 be resolved on our grid. 

— Specify the temperature 

{-^core 
if Pi_i > pHe 
Tcore + |(Tbasc " T.orc) [l + tauh (M^)] if i [l + tauh (^^^^f^)] < 0.999 
max{T(sHe, Pi, ^i), ^cutoff} otherwise 

(2) 

here, Tbasc is the desired temperature at the base of the He layer, and suc is the 
specific entropy of the base, sue = s(pbase, ^base, -^basc)- Finally, Tcutofr is the lowest 
temperature allowed in the outer envelope. 

In both of these profiles, a tanh profile was used at the base of the layer. The 
true functional form of the transition is not known. Adopting a tanh gives it 
a smoothness that is desirable for hydrodynamics codes while still keeping the 
transition narrow. 

— Compute the pressure using the equation of hydrostatic equilibrium (HSE). We 
difference the HSE equation as: 

Pf^"" = Pi~i + l^Pi + Pi^i)9i-i/2 (3) 

Here gi-1/2 is the gravitational acceleration at the lower edge of the zone, com- 
puted as gi^i/2 = —GMcnci/ x1_^i2, where Mend is the mass integrated up to that 
edge. Given our guess for the density in the zone, p^, the equation of state will re- 
turn a pressure pf°^ = p(pi, Tj, Xi). We define a function, F, as F = pf^^ — pf*-*^ 
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and use a Newton-Raphson iteration to find the pi needed to zero F. In the case 
where we are constraining the profile to be isentropic (with constant entropy shc) 
we begin with a guess for Tj, and zero an additional function, G = sue — sf*-*^, 
where sf^^ = s{pi,Ti, Xi). In this case, zeroing both F and G yield pi and Tj. 

— After all zones are computed, compute Mwd by integrating up all the carbon (and 
oxygen if present), and Mne by integrating up all the helium. We then correct pcore 
and phc using a secant method, recompute Pcore, and iterate the above procedure 
until the model is converged. 



For the current simulations, we use two slightly different models. For both, we set Xcore 
to be pure carbon and to be pure helium. (The reason we leave oxygen out of the white 
dwarf for these initial models is to allow us to easily use the oxygen generated by the reactions 
as a tracer of the nucleosynthesis). We pick Mwd = 1 M© (here, Mwd is the integrated mass 
of the carbon only), Mrc = 0.05 Mq, Tcorc = 10^ K, Tbasc = 2 x 10^ K, and 5 = 5 x 10^ cm. 
This value of 6 ensures that the composition gradients are somewhat smoothed. The only 
difference between the models is the choice of cutoff temperature above the convective zone. 
Our cool model has Tcutoff = 5 x 10^ K, and our hotter model has Tcutos = 7.5 x 10^ K. 
Figure [1] shows the profile for both models. The cool model is indicated by the solid lines 
and the hot model is the dashed lines. When creating these models, we pick a zone width, 
Ar, to be 1/5*^ of the Cartesian zone width. Ax. 

It is interesting to look at the behavior of the sound speed in the initial model, shown in 
the lower panel of Figure [H For very low Mach number flows, the time step for a compressible 
algorithm (assuming uniform zoning) will be set where the sound speed is highest — the center 
of the star. For the low Mach number algorithm, it will be set where the velocities are the 
highest — presumably in the convectively unstable region. But the highest Mach number in 
our simulations is not necessarily there as well — it is likely to be at the very edge of the 
star, where the velocity will rise due to the density gradient, but the sound speed is small 
relative to the sound speed at the core. For this reason, we may realize a moderately large 
Mach number just outside the star (~ 0.2-0.3), but still be able to take a time step an 
order-of-magnitude larger than a compressible code, because the peak Mach number is not 
where the sound speed peaks. We note that the upturn in sound speed at the largest radii 
shown in Figure [1] is not mapped into our computational domain because of the use of a low 
density cutoff. This upturn arises because of the dominance of the radiation pressure in the 
general equation of state, a behavior that is unphysical outside the star. 
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2.3. Simulation Parameters 



W e use the same s et of c utoff densities and sponging tecliniques described in lZingale et al. 



( 120111 ) : lNonaka et al.l (120 12[ ). These parameters are designed to ehminate large velocities that 
arise due to the steep density gradient at the edge of the star, outside of our region of interest. 

A low density cutoff, Pcutofr, is the minimum density used in the initial model — outside 
of that radius in the star, rcutoff = ''"(Pcutoff), we hold the density constant. The material 
at r > Tcutoff does not contribute to the computation of the gravitational acceleration. To 
determine a reasonable value of Pcutoff? we perform a one- dime nsional base state exp ansion 
test, comparing the results to the compressible code CASTRO (lAlmgren et al.ll2010bl ). Fig- 
ure [2] shows the results of placing our initial model onto a one-dimensional spherical grid 
and heating it for 25 s with a heating term: 



H = AheatX(^He)e 



-{r-rhe3,t)'/cr\ 



2 

heat 



(4) 



We choose Ahcat = 10 erg g~ s" , rhcat = 4.2 x 10 cm, and aheat = 10 cm. These values 
were chosen to concentrate the energy release at the base of the He layer and to ensure 
that we see a significant response to the hydrostatic s tructure. The base state expansion 
algorithm used in MAESTRO is described in detail in iNonaka et al.l (120 lOf ) . We run with 
three different choices of PcutoS : 5 x 10"^, 10^, and 1.5 x 10^ g cm~^. The figure shows excellent 
agreement between the fully compressible (CASTRO) results and MAESTRO's base state 
expansion. We note that this is a rather severe test, and the amount of expansion seen 
here is greater than what we expect in our three-dimensional simulations. There is a slight 
departure from the compressible solution at the base of the helium layer in the temperature 
field for the largest choice of PcutoS- Bas ed on these results, we choose Pnitnff = 10 '*' g cm~^ 



for the simulations presented here. As in lZingale et al.l (1201 if ): iNonaka et al.l (120 12[ ) we have 
a cutoff to the buoyancy term in the momentum equation which we set as 2pcutoff- It is 
important to note that this one- dimensional test does not have any transport of the energy, 
it simply expands the hydrostatic structure in response to the heating. As a result, the large 
increase in temperature seen here will not appear in the actual three-dimensional simulation 
(since convection would redistribute the heat). 

Outside of the star, we want to damp away large velocities, as this region is not really 
part of the simulation space and we do not want velocities here to build up and influence our 
time step. In MAESTRO, this is accomplished through the use of a sponge, which appears 



in IZingale et al.l (120091 ) and as refined in IZingale et al.l ( 120111 ): iNonaka et al.l (120121 ). In 
particular, we define a density at which to center the sponge, p^d, and a multiplicative factor, 
/sp, to mark the start of the sponge. This means that the sponging turns on (gradually) once 
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the model drops below a density of /spPmd- We chose the sponge parameters to make the 
sponge turn on at the top of the convective layer. Finally, we modify the constraint equation 
to act li ke the anelastic const raint to suppress high velocities at the outer boundary of the 
star (see lAlmgren et al.ll2008[ ). This occurs once the density drops below Paneiastic- Tabled] 



summarizes the cutoff parameters for the two different initial models. These parameters are 
illustrated in Figure [1] as the vertical lines, again with the solid corresponding to the cooler 
cutoff model and the dashed to the hotter cutoff model. 



2.4. Grid Structure 



Adaptive mesh refinement is used, with the refinement tagging on zones that have 
X(^He) > 0.01 and p > Pcutoff- Additionally, we always refine the very center of the star (the 
coordinate origin) du e to the design of t he averaging algorithm from the Cartesian grid to 
the radial base state (INonaka et al.ll2010l ). Simulations are run with one level of refinement 
with a factor of two increase in resolution on top of the coarse grid — Figure |3] shows a 
representative grid structure. The work is parallelized by distributing grids to nodes that 
communicate with each other using MPI and OpenMP to spawn threads within nodes to 
perform floating point work on the data. Unless otherwise noted, all simulations use a base 
grid of 256^ with an additional level of refinement on the helium layer. The computational 
domain is a cube with a side of 7.5 x 10^ cm, with the center of the star placed at the origin. 
This gives a resolution of Ax'^'^° = 14.6 km in the convective region. The one- dimensional 
radial grid for the base state uses a resolution Ar = to improve the performance of 

the mapping between the Cartesian and radial grids (see lNonaka et al.ll2010l ). The boundary 
conditions are reflecting on the symmetry faces of the domain (lower x, y, and z) and outflow 
(zero- gradient) on the other faces. 



2.5. Velocity Field Initialization 



We d efine the initi a l velo city field to be a perturbation with similar form to that de- 
scribed in lZingale et al.l (120091 1. For a zone with coordinates {x, y,z), a set of Fourier modes 
is defined as: 



*-^Z,m,n — COS 



Lm,n 



iv) 

m,n 



COS 



27imy 



a 



Lm,n 



c, 



(^) 

l,in,n 



COS 



2nnz 



+ 



a 



l,m,n j ' 



(5) 
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Lm,n I ' Lm,n 



(y) 



sm 



dy) 



a 



Lm,n I ' Lm,n 



sm 



f 2^ + 



l,m,n j ' 

(6) 



where a is the perturbation scale and 4>fmn^ randomly generated phases that lie between 
and In. The velocity perturbation in the zone is then set as: 

3 3 3 



l=\ m=l n=l 
3 3 3 

EEE 

1=1 m=l n=l 
3 3 3 

XI XI XI AT, 



1 

'n, 

' l,m,n 
1 



- n,m,nl 1 1'^ i^rn,n'^l,m,n^ l,m,n yl,m,nlf^i^rn,n'^l,m,n'-'l,Tn,n 



jq{x) My) Mz) _ r,r^^> '^^''> 



{^) My) c(^) 



(7b) 

i=l m=l n=l 

where the amplitudes a/,m,n, (^i,m,n, and 7/,m,n are randomly generated to lie between —1 
and 1, and Ni^rn,n = a//^ + + n^. 

Finally, we confine the perturbation to lie in the convective region as: 
Au' I r^"':X" — r — a\ I r — r 



u 



4 

Av' 
Aw' 



w 



1 + tanh 
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r — 
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1 + tanh 
1 + tanh 
1 + tanh 
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-d 



d 



pert 



-d 



d 

inner J 

pert " 



d 



(8a) 
(8b) 
(8c) 



Here, the first tanh factor cuts the perturbation off at the outer edge of the star and the 
second tanh factor cuts the perturbation off at the base of the convective boundary. These 
transitions are characterized by a thickness d. A is the amplitude of the perturbation. 
For all the simulations presented here, we choose A = 10^ cm s~^, a = 5 x 10^ cm, and 
d = 10^ cm. The inner extent of the perturbation, r™^'^^ is set to be the radius where 
X('*He) > 0.9 (moving from the center of the star to the edge). The outer extent is set as 
'^perr = ['^perr + '"(/spPmd)] /2, wherc ^(/spPmd) is the radius where the sponge just begins 
to turn on. This confines the velocity perturbation to the lower half of the convectively 
unstable layer. 



2.6. Hydrodynamic Integration Strategy 



The construction of advective fiuxes on the faces of the computational zones requires 
the prediction of values of the fiuid state from the cell centers to the faces at intermediate 
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times. The method described in iNonaka et al.l ( 120101 ) predicts the mass fractions, X^, and 
the perturbational density, p' (the difference between the zone's density and the base state 
density at that radius), to the interface and algebraically combines these with the base state 
density to compute (pXfe) on the edges. Here, we instead predict the full density, p, and 
and define the edge state as the product of these. Numerical experiments show that this 
variant is more robust at steep composition gradients (such as the gradient at the base of 
our helium layer), in that unperturbed gradients become les s smoothed. This is similar to 
the treatment of the density in our original implementation ( lAlmgren et al.ll2006l ). 



3. Results 



Unlike our previous simulations of interior convection ( IZingale et al.ll2011t iNonaka et al. 



2OI2I ). in this problem convection takes place on the surface of the star. Here we will assess 
the applicability of our simulation methodology to modeling the convection in the helium 
layer on the surface of the white dwarf. We discuss the general trends of the simulations 
and explore the effect of various simulations parameters. We defer a detailed discussion of 
the scientific implications and explorations of other white dwarf masses to a later paper. 

Our reference calculation is an octant simulation with the hot-TcutoS initial model. We 
will compare results from this reference calculation to a simulation with: 



disabled burning (Section 
the full star rather than an octant (Section 13. 3p . 
twice the spatial resolution (Section 13. 4p . 

the cool-Tcutoff initial model rather than the hot-Tcutofr model (Section l3.5p . 

a modified treatment of the region beyond the convective layer, as controlled by the 
cutoff densities and sponging parameters (Section l3.6p . 



3.1. General Trends 

We begin by looking at the qualitative behavior of the convection for our reference 
calculation. Each time step we store the location of the hottest zone in the convective region 
(only considering cells with p > fspPmd)- We also store the maximum Mach number in the 
entire computational domain. 
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Figure H] shows a time sequence for the reference calculation, with both the ^^O abun- 
dance and radial velocity (red = outflow; blue = inflow) visualized. The development of the 
convection is clearly seen. By 80 s, the ^^O synthesized at the base of the layer is distributed 
throughout the entire convective zone, and a clear top to the shell is seen. The radial velocity 
shows that the convection is divided into cells, with outflow regions surrounded by inflowing 
regions. Looking at the 80 s plot, there are approximately 6 such cells along the edge of the 
domain. The circumference of this edge is ~ (7r/2)5 x 10^ cm or 7.9 x 10^ cm. This tells 
us that a convective cell has a diameter of approximately 1.3 x 10® cm. This is very close 
to the thickness of the convective layer (~ 10® cm as seen in Figured]). We note that the 
thickness of the convective layer is a function of the initial model, and we will explore other 
initial models in subsequent papers. 

Figure shows the peak temperature and peak Mach number as a function of time 
for this run. It is interesting to see that there are a number of failed ignitions toward the 
end — our suspicion is that the buoyant hotspots rise and cool via expansion faster than the 
rate at which energy is injected through reactions. This is supported by the relic plumes of 
^^O seen throughout Figure HI A caveat to this behavior is that our network only continues 
to the production of ^^O. A more extensive reaction network might release enough energy 
for the first hot spot to fully ignite. 

The Mach number panel in Figure [5] shows that the Mach number stays below 0.3 for 
the bulk of the simulation. It is important to note that the efficiency metric for the low Mach 
number algorithm, the increase in time step over a compressible code, is not simply 1/M here. 
For a compressible simulation, the higher sound speed at the center of the star dominates the 
fluid velocities realized at the edge of the star. For these calculations, we evolve to 100 s in 
~ 2100 time steps, giving an average AtiowMach ~ 0.05. The peak sound speed in these models 
is Cs = 4.7 X 10® cm s~^. Taking the coarse grid spacing. Ax = 2.9 x 10^ cm, the compressible 
time step for these models would be At = Ax/cs = 6.2 x 10^'^ s — a factor of 8 smaller. Note 
that this doesn't factor any of the fluid motions themselves into the compressible time step, 
which would only further reduce it. We also note that the early evolution begins with a 
much smaller Mach number, so the efficiency is greatest earlier in the evolution. Finally, 
once a hotspot ignites, the Mach number rises beyond the range of validity of the low Mach 
number approximation, and we would need to transition this problem to a compressible code 
to continue. Initial w ork on transitioning MA ESTRO calculations to the compressible code 



CASTRO is shown in lAlmgren et all fl2010ar ). 



Figure shows a side view of the convective layer for the reference calculation. We see 
that the convective layer is well bounded. The grey contour marks the density at which the 
sponge just begins to turn on — as designed, this is at the upper boundary of the convection. 
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The orange and green contours show the regions where the Mach number is highest. We see 
that these regions are very small points, scattered throughout the layer, mostly near the top 
of the convective zone. 

Figure [7] shows the location of the hottest point in the model (the hotspot) at each time 
step, colored by the logarithm of the temperature. We note that the time step changes over 
the course of the simulation (generally speaking, smaller time steps toward the end of the 
simulation). The points where ignition developed (seen, for example, as the ^^O plumes in 
the final time snapshot in Figure H]) are shown in deep red. We notice that there seems to be 
a slight bias of ignition along the edges of the octant, where the flow can be confined by the 
geometry. We remind the reader that all simulations are done in Cartesian coordinates — the 
meridians and parallels drawn on the sphere are for guidance only. 

At a single snapshot, it is interesting to see how many hot plumes exist. Figure [S] shows 
the distribution of the hottest zones within the p — T plane at two different times. The 
color coding of the dots indicate spatial location within the finest level of refinement on 
our grid — in particular the normalized x,y, z coordinates of each zone are translated into 
r,g,b colors for each data point. Groups of points of a similar color indicate simulation 
zones in close spatial proximity. The left plot, at t = 98 s, shows in purple a single plume 
dominating the flow. There are a few smaller regions beginning to heat, the largest of which 
is the clump of green points. Only four seconds later [t = 102 s), as seen in the right plot, 
the previously dominant plume (purple of the left plot) has cooled and redistributed its 
heat while several other distinct plumes have formed. This indicates that there are possibly 
several ignition regions, or at least several regions that are almost to runaway conditions 
once the dominant hot spot (blue points in this case) ignites. The extreme temperature 
sensitivity of the reaction rates makes it difficult to determine if these other "failed" ignition 
points would actually ignite. 

Finally, we can look at how much the atmosphere expanded over the course of the 
simulation. Figure [9] shows the base state density at the start of the reference hot-Tcutoff 
simulation and after 100 s of evolution, as well as the temperature averaged in a shell of 
constant radius at the start and end of the simulation. We see that a bit of expansion 
has taken place, smaller in magnitude but qualitatively the same as in our one-dimensional 
test. This supports the statement that we need to use the base state expansion in the 
hydrodynamics to accurately model the flow. The temperature plot shows that the heat 
generated by reactions has been transported throughout the convective layer. We note that 
there is heating below the convective layer, which may be due to the adjustment from the 
expansion of the base state on the adaptive grid. The magnitude of this heating is small and 
we do not expect it to affect the results. This is not seen in the case where we do not burn 



(next section). 



3.2. Disabled Burning 

To demonstrate that the convective velocity field we see is driven by the reactions 



and not due to the initial model's temperature profile (IMocak et al.l 120101 ) or discretization 
error, we run a test with the initial velocity perturbations, but burning disabled. The radial 
velocity field in this case is shown in Figure [TOl We note that the range used in the contours 
is smaller than in Figure S] to bring out the detail. We see no suggestion of the convective 
pattern that dominates in the burning calculations. Figure E] shows the temperature and 
Mach number over the course of the simulations for the reference calculation with burning 
disabled. We see that the peak temperature stays right around the starting value of 2 x 10* K, 
as expected. The peak Mach number hits a plateau of 0.05, driven by the artificial buoyancy 
introduced by the mapping error between the one-dimensional hydrostatic base state and the 
three-dimensional Cartesian state. This value does not seem to grow further. It is important 
to note that a compressible code would also see a velocity field generated in this test, again 
driven by the inability to exactly satisfy hydrostatic equilibrium numerically. Taken together, 
these figures indicate that the convective behavior described above arises due to the energy 
release from the reactions. 



3.3. Full Star Simulation 

We run a single full-star calculation (using the hot-Tcutofr model) to assess the infiuence of 
the octant geometry on the general results described above. The resolution is the same, with 
the base grid now twice as large in each coordinate direction. The convective field is shown 
in Figure [TTl and we see the same overall structure that appears in the octant simulations. 
The time-dependent peak temperature and Mach number (Figure [T2|) also agree well with 
the octant case. Figure [13] shows the hotspot location over time for this calculation. We see 
a uniform distribution of points over the sphere. 



3.4. High-Resolution 

To understand the robustness of the convective features to resolution, we run a single 
case (hot-TcutofT model) at twice the resolution. This is accomplished by doubling the number 
of zones in each direction on the coarse grid. Figure [14] shows the convective field for this 
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simulation, in which the overall structure agrees with that of the standard-resolution simu- 
lations. The time- dependent peak temperature and Mach number for this run (Figure [T2l) 
again show excellent agreement with the standard-resolution runs. This gives us confidence 
that we are sufficiently numerically converged. 



3.5. Cooler Initial Model 

The corresponding images for the cool-Tcutofr simulation near the point of ignition are 
shown in Figure [TH There is a strong qualitative similarity with the hot-Tcutofr run, indicating 
that the structure of the convection is insensitive to the details of the top of the convective 
zone. We also see the expected behavior that the hot-Tcutoff model has a lower peak M 
than the cool-Tcutoff model because of the higher sound speed at the edge of the star. This 
suggests that the peak Mach number occurs at the edge of the star, outside of the convective 
region. 



3.6. Effect of Cutoff Densities and Sponging 

In the simulations above, we chose the sponging parameters such that the sponging 
begins just at the top of the convective layer. Here we explore the sensitivity of that choice 
by turning the sponge on lower — we now set = Panciastic = 10^ g cm~'^. Additionally, we 
decr ease the timescale oy er which the sponging acts by changing /t = 10 s~^ to k = 100 s~^ 



see 



Almgren et al.l 120081 ). 



Figure [16] shows the structure of the convective field. We notice that the radial extent 
appears slightly diminished compared to the previous simulations, owing to the more ag- 
gressive sponging. However, the overall pattern of convective cells appears consistent with 
the other cases. The trends of peak temperature and Mach number are shown in Figure [TH 
and exhibit a slightly lower peak Mach number as the reference simulation due to the more 
aggressive sponging. These comparisons show that the convective behavior is not strongly 
dependent on how we treat the top of the convective layer. 



4. Conclusions and Discussion 

The main goal of the present paper is to serve as a proof-of-concept that the low Mach 
number methodology can be applied to shell burning on the surface of white dwarfs. This is 
the first application of MAESTRO where we have had off-center heating with an expanding. 
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self-gravitating hydrostatic state. We demonstrated that efficient three-dimensional models 
of the convective flow leading up to the ignition of a burning front in a helium layer on 
the surface of a white dwarf are possible. We explored the sensitivity of our results to a 
variety of factors and showed that the convective features realized are robust. The octant 
calculations are inexpensive to run (requiring only about 40,000 CPU hours on the OLCF 
jaguar machine, using 128 MPI tasks with 8 OpenMP threads per task). This suggests that 
a parameter study of progenitor models is feasible — this will be the focus of a follow-on 
paper. There is a wide variety of potential models — varying white dwarf masses and helium 
envelopes. The conditions at the base of the helium layer will vary across these different 
models, so some models may be more amenable to our methodology than others. This needs 
to be explored on a case-by-case basis. For instance, we expect that low mass shells and 
lower mass white dwarfs would have slower dynamics. 

Future work will focus on better understanding the conditions in the helium layer lead- 
ing up to ignition. Open scientific questions that we wish to understand are whether the 
ignition occurs in a manner that is amenable to a detonation — this is a key requirement for 
the sub-Chandrasekhar explosion models. Also of interest is whether ignition can arise in 
multiple disjoint points on the surface of the white dwarf. Studying this will require either 
enhancements to the low Mach number model (i.e. the addition of long wavelength acoustics) 
or feeding MAESTRO models into a compressible hydrodynamics code. Finally, the models 
presented here started rather late in the evolution, but it is easy to start with cooler models 
to see more of the ramp up to ignition. This earlier evolution will take place at lower Mach 
numbers. 



Videos of the reference calculation are available at: http://youtu.be/boHVbcfazvw 



and http://youtu.be/37WqQ0Km0p4 We thank Frank Timmes for making his equation of 
state routines publicly available and for helpful discussions on the thermodynamics. The 
work at Stony Brook was supported by a DOE/Office of Nuclear Physics grant No. DE- 
FG02-06ER41448 to Stony Brook. The work at LBNL was supported by the Applied Math- 
ematics Program of the DOE Office of Advance Scientific Computing Research under U.S. 
Department of Energy under contract No. DE-AC02-05CH11231. 

An award of computer time was provided by the Innovative and Novel Computational 
Impact on Theory and Experiment (INCITE) program. This research used resources of the 
Oak Ridge Leadership Computing Facility located in the Oak Ridge National Laboratory, 
which is supported by the Office of Science of the Department of Energy under Contract 
DE-AC05-00OR22725. Visualizations were performed using Visit and matplotlib. 
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Fig. 1. — The initial model. The solid lines are for the cooler cutoff model (Tcutofi = 5 x 10^ K) 
and the dashed lines are for the hotter cutoff model (Tcutos = 7.5 x 10^ K). The vertical lines 
indicate the start of the sponge (leftmost), the anelastic cutoff (center), and the base cutoff 
density (rightmost), again with solid for the cooler cutoff model and dashed for the hotter 
cutoff model. 
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Fig. 2. — A test of the base state expansion algorithm. The initial model is shown as the 
dotted line (for the choice Pcutofr = 10"^ g cm~^). The compressible solution by the CASTRO 
code is shown as the black line, and three MAESTRO solutions are shown, with varying 
Pcutoff- We see very good agreement across the board between MAESTRO and CASTRO, 
with the two lowest cutoff densities showing the best match. 
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Table 1. Cutoff parameter values 



parameter 


cool-TcutofF model 


hot-Tcutoff 


mo 


del 


Pcutoff 


10^ g cm~^ 


10^ g cm" 


-3 




Pmd 


3 X 10^ g cm~^ 


6 X 10^ g 


cm" 


-3 


fsp 


2.0 


2.0 






Panelastic 


3 X 10^ g cm-3 


6 X 10^ g 


cm" 


-3 




Fig. 3. — A representative 2-level grid used for the simulations. 
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Fig. 4. — Time sequence of the ^^O abundance (left) and radial velocity (right) for the hot- 
^cutoff reference simulation shown at 20, 40, 60, 80, and 100 s. The tick marks are spaced 
10^ cm apart. The ^^O sequence highlights the extent of the convective region well — by 80 s, 
there is a well-defined upper surface to the convective zone. The radial velocity (CGS units); 
red=outflow, blue=inflow) shows a characteristic "granulation" pattern, making the distinct 
convective cells stand out. 
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Fig. HI — Time sequence continued. 
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Fig. 5. — Peak temperature and Mach number as a function of time for our reference calcu- 
lation (hot-Tcutoff) and the same calculation with burning disabled. 
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Fig. 6. — Side view of the convective field of the hot-Tcutofr reference run at 95 s. All contours 
are in CGS units. The tick marks are 10^ cm apart. The grey density contour marks the 
surface where the sponging just begins to turn on — we see that the convection is entirely 
below this surface. The orange and green contours show the locations where the Mach 
number is highest. 
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Fig. 7. — Location of the hottest point at each time step in the reference simulation, colored 
according to temperature. The red boundary shows the bounds of the octant modeled. 
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le9 ^ ^ = 98 S ^ ^ ^ j^ le9 ^ t = 102 5 




P(gcm-'') p(gcm-') 

Fig. 8. — Distribution in the p — T plan of the hottest zones in the finest level; the left 
plot shows a snapshot at t = 98 s and the right plot shows a snapshot at t = 102 s. Colors 
represent a proxy for spatial location. Within four seconds, the single dominant plume in 
the left plot is washed away and replaced by several distinct hot regions of which the group 
of blue zones reach runaway conditions. 
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Fig. 9. — Base state density and lateral average of the temperature at the start and end of 
the reference hot-Tcutoff simulation. 
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Fig. 10. — Radial velocity (CGS units) for the "no burning" test model after 197 s. Note: 
the range on the contours is smaller here than in Figure H] to bring out detail. The tick 
marks are 10^ cm apart. We see no evidence of the convective field structure. Instead, the 
random velocities here are driven by discretization error. 
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Fig. 11. — ^^O abundance (left) and radial velocity (right; CGS units) for the fuUstar hot- 
^cutoff niodel at 96 s — ^just at the point where we are igniting. The tick marks are 10^ cm 
apart. The overall convective structure compares well to the octant simulations. 
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Fig. 12. — Peak temperature and Mach number as a function of time for the various sup- 
porting calculations: full star, high-resolution, cool-Tcutofr model, and stronger sponging. 
The reference hot-Tcutofr calculation is also shown. We see that the trends for all quantities 
are the same for all runs, indicating that our simulation is converged in resolution, and the 
octant is a reasonable model for the energetic evolution, and that the treatment at the outer 
radius of the convective region is not critically important. 
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Fig. 13. — Location of the hottest point at each time step for the full-star, hot-Tcutofr calcu- 
lation. 
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Fig. 14. — ^^O abundance (left) and radial velocity (right; CGS units) for the high-resolution 
hot-Tcutoff near the point of ignition (98 s). The tick marks are 10^ cm apart. The overall 
structure compares well to the standard-resolution hot-Tcutofr simulation. 
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Fig. 15. — ^^O abundance (left) and radial velocity (right; CGS units) for the cool-Tcutoff 
model at 100 s. The tick marks are 10*^ cm apart. The overall structure compares well to 
the hot-Tcutoff simulation. 
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Fig. 16. — ^^O abundance (left) and radial velocity (right) for the hot-Tcutofr run with more 
aggressive sponging at the top of the convective layer. Shown at 100 s. The tick marks are 
10^ cm apart. Again, we see good agreement with the other cases. 



